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Abstract 

Introduction. The task of analyzing the stability of plates and shells under creep conditions is critical for structural 
elements made of materials with the property of aging, which are under the action of long-term loads, since the loss of 
stability can occur abruptly and long before the exhaustion of the strength resource of the material. Currently, the issues 
of joint consideration of geometric nonlinearity and creep in the problems of buckling plates remain poorly studied, 
existing software systems do not provide such calculations. The objective of this work is to develop an algorithm for 
calculating the stability of rectangular plates with initial deflection, which are subjected to loads in the middle plane, 
taking into account geometric nonlinearity and creep. 

Materials and Methods. When obtaining the resolving equations, the geometric and static equations of the theory of 
flexible elastic plates were taken as the basis. Physical equations were derived from the assumption that total strains were 
equal to the sum of elastic strains and creep deformations. Finally, the problem was reduced to a system of two differential 
equations, in which the desired functions were the stress and deflection functions. The resulting system of equations was 
solved numerically using the finite-difference method in combination with the method of successive approximations and 
the Euler method. As the boundary conditions for the stress function, the frame analogy was used, as in the case of a plane 
problem of elasticity theory. 

Results. The solution to the problem for a plate compressed in one direction by a uniformly distributed load has been 
presented. The nature of the growth of displacements at different load rates and initial deflection was studied. It has been 
established that when the vertical displacements reach values comparable to the thickness of the plate, their growth rate 
begins to decay even at a load greater than the long-term critical one. 

Discussion and Conclusion. The results of stability analysis using the developed algorithm show that the growth of plate 
deflection under the considered boundary conditions is limited, stability loss is not observed at any load values not 
exceeding the instantaneous critical one. This indicates the possibility of long-term safe operation of such structures with 


a load less than instant critical one. 
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AHHOTalna 


Beedenue. 3ajaua ananv3a yCTOMYMBOCTH MIaCTHH HW OOOMOUeK B YCOBHAX MOU3YYeCTH aKTYAaJIbHa [JIA JIEMCHTOB 
KOHCTpyKIUHii H3 MaTepHasioB, OOaqaloulMx CBOMCTBOM CTapeHHA, HaxXOJAMWIMXCA HO WeviCTBHeM JJIMTeJIbHBIX Harpy30K, 
HOCKOJIbKY HOTepx YCTOMYMBOCTH MOXKET MpOHCXOAMTb pe3KO M 3aONrO WO UCueplaHHaA MpOuHOCTHOTO pecypca 
MatepHana. Bompocbl coBMecTHOrO y4eTa TeOMeTPHYeCKOM HeMHeEMHOCTH MU MOU3y4ecTH B 3aayax BbINyYMBAaHHA 
TUIACTHH B HAaCTOAINee BPEMA OCTAIOTCA CIaOO H3YYeHHBIMH, CYWIECTBYIOMMe IPOrpaMMHBIe KOMIJICKCHI He MO3BOJIAIOT 
BBINOJHHTb Tako pacuér. Llenb1o HacTosMel paOoTEI BLICTyMaeT pa3spaOoTKa ajIropuTMa pacueTa Ha yCTOMUYMBOCTb 
TIPAMOYTOJIBHBIX TIACTHHOK C HavaIbHOM MOTHOBIO, HCIIbITHIBAIOMIMX elicTBHe Harpy30K B CpeJ{MHHOM MIOCKOCTH C 
yu4eTOM TeOMeTPHYeCKOM HeIMHeHHOCTH HW MOU3y4ecTH. 

Mamepuanei u memoovi. pu nonyyenuu pa3pelmaroulMx ypaBHeHHii B OCHOBY MOJO%#KeHbI TeoMeTpHyeckHe 
cTaTHueckve ypaBHeHHaA TeOpuu rHOKHX yUpyrHux WiacTHH. Du3st4eckHe ypaBHeHHA BbIBOATCA 13 MpeANOOKeHHA, 4TO 
HOMHbIe WepopMalHu paBHbI CyMMe yupyrux WedopMalui uv epopMalui Nom3yyecTH. OkoHYaTeIbHO 3ayjaya Obuia 
cBeyeHa K CHCTeMe 13 AByx AuddepeHMasbHbIX ypaBHeHHH, B KOTOPbIX B KaYeCTBE HCKOMBIX (PyHKI[MM BbICTYHaIOT 
(YHKUMA HallpwKeHH u npornOa. PemieHue NouyYeHHOM CHCTeMBI ypaBHeHHM BBITIOIHAIOCb YMCIICHHO C MOMOLIbIO 
MeTOJja KOHCYHBIX pa3HOCTell B COUeTAHHH C MeTOZOM MOCIeAOBaTeIbHBIX IpHOMMKeHHM HU MeTOZOM Diimepa. B 
Ka4eCTBe TpaHHYHbIX YCNOBU Id PYHKUMM HallpsxKeHHM MCHONL3yeTcA PpaMHad aHaJIOrMA, KaK B CJIydae TMIOCKOM 
3aja4H TEOpHH yIpyrocru. 

Pezyismamot ucciedosanua. B pamkax TlocTaBieHHO Wem paspadotaH aIropHTM pacueta H pescTaBeHo pelleHue 
3aaqH JIA TWIaCTHHbI, CAKMMACMOM B OJHOM HallpaBJIeHHH paBHOMepHO paciipeyeneHHo Harpy3Kkou. UccneqoBan 
XapaKTep pocta MepeMeleHH pH pa3IM4HOM BeIMYMHe Harpy3KH HW HayabHOH WorMOu. YcTaHoBeHoO, 4TO pu 
JOCTWKEHHM BeEPTHKAIbHbIMH TepeMeLICHHAMH BCIMYHH, COM3MEPHMBIX C TOJIMIMHOM MWIaCTHHKU, CKOPOCTb HX pocta 
HauyHHaeT 3aTyXaTb Jake IPH Harpy3Ke OobUe JIMTebHOM KpUTHYECKON. 

O6cystcdenue u 3akiiouenue. Pe3ynbTaTbI aHaliu3a yYCTOMYMBOCTH C MCMOb30BaHHeM pa3paOoTaHHorO aJiIroputma 
MOKa3bIBaloT, YTO pocT MporHOa WiacCTHHbl pH PaCCMOTPeHHBIX TpaHW4HbIX YCJIOBHAX OTpaHHyeH, WoTeps 
yCTOM4MBOCTH He HaOOaeTCA Mp JNOObIX 3HAYCHHAX Harpy3KH, He MIpeBOCXOMAWHX MFHOBCHHY!IO KPHTHYECKyIO. ITO 
TOBOPHT 0 BOSMOXKHOCTH AJIMTeIbHOM Oe30nacHol SKCIyaTalMu TaKHX KOHCTpyKUMi pu Harpy3ke MeHee MrHOBeHHOM 
KpHTH4eCKOH. 


Kor0ueBble CJIOBa: yCTOH4MBOCTB, TOJI3YU4CCTh, Wi1acTHHa, rCOMeTpH4eCcKaAa HeJIMHeCMHOCTS, dbu3sn4eckaa HeJIMHCHHOCTS, 
HadaJIBHbIle HCCOBCpIICHCTBa, MCTO, KOHCUHBIX pa3HocTeii 


BuaarojapHocrnu: aBTOPbI BbIpaxKaloT OmarolapHOcTE peqakiHv HU peleH3eHTaM 3a BHHMAaTeCJIBHOe OTHOMICHHE K CTaTbe 
HW yKa3aHHble 3aMedaHhA, KOTOPbIe MO3BOJIMJIM HOBbICHTb Ce KadecTBO. 


Aaa wntTupopanna. A3piepn C.b., UenypHenxo A.C. Bsmyanpanve pAMOYrOJIbHbIX TWiaCTHH TIpH HesMHeMHOH 
mom3yuectu. Advanced Engineering Research (Rostov-on-Don). 2023;23(3):257-268. https://doi.org/10.23947/2687- 
1653-2023-23-3-257-268 


Introduction. Much attention is paid to the stability analysis of thin-walled structures in the form of plates and shells, 
since such structures are widely used in construction and other branches of technology [1-3]. One of the challenges in 
the field of calculating plates and shells is the analysis of their stress-strain state under creep conditions, which is 
confirmed by a significant number of works published recently on this problem in domestic and foreign sources. Thus, 
in [4-8], the issues of buckling under creep of composite thin-walled structures were investigated. In [9], the problem of 
stability of functionally gradient plates was considered, taking into account the dependence of material properties on 
temperature. In [10], stochastic analysis methods were applied to the problem of buckling composite plates. In [11-17], 
the issues of stability of viscoelastic plates and shells under the influence of dynamic and tracking loads were discussed, 
and [18] dealt with plates of medium thickness, taking into account the dependence of material properties on time. 
Mathematical difficulties arising in solving these problems led to the fact that numerous researchers limited themselves 
to linear laws of viscoelastic deformation or considered the case of steady-state creep. The finite element method opens 
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up great possibilities in solving the problems of calculating plates and shells taking into account creep. However, modern 
computational complexes, such as ANSYS, Abaqus, LIRA, etc., contain a limited set of rheological models applicable to 
specific materials in a fixed range of stresses and temperatures. There is a need for alternative calculation methods suitable 
for arbitrary laws of viscoelastic deformation, including nonlinear ones. 

This work was aimed at constructing a system of resolving equations for the problem of buckling rectangular plates 
with nonlinear viscoelastic properties under the action of forces in the middle plane, taking into account large 
displacements, as well as an algorithm for its solution. Note that the problem of the stability of structural elements, taking 
into account creep, cannot be solved as a problem of pure stability. Its solution requires the presence of disturbances in 
the form of initial irregularities. Generally, the initial imperfections are given in the form of the initial loss or eccentricities 
of the application of loads. 


Materials and Methods. Let us consider the calculation method on the example of a plate with an initial deflection 
Wo(x,v), compressed by a distributed load p [kN/m] in the x-axis direction and having a hinge support along the 
contour (Fig. 1). 


Fig. 1. Computational scheme 


In the case under consideration, in the presence of creep, if compared to the theory of elastic flexible plates, the 
difference manifests itself only in the form of physical equations. Total deformations can be presented as the sum of 
deformations of the middle plane (passing into the surface) and bending deformations caused by changes in the curvature 


of the middle surface: 
02w 02w 02w 


Vx =Y°-22 ; (1) 


where €, and ¢, — total linear deformations; y,,— total angular deformations; ¢° and ¢° — linear deformations of the 


middle surface; y° — angular deformations of the middle surface. 


For deformations of the middle surface, the equation of continuity of deformations can be written [19]: 


076° eo Ory? _ | (w+ wy) > @ (w+ Wo ) CO? (w+ Wy ) Ow, . O2W, 07 Wy Q) 
Oy? Ox? Oxy Oxdy Ox? oy? Oxdy Ox? dy? | 
For materials with viscoelastic properties, total deformations can be represented as: 
l ee Sige au 
£.= ml ove, ) +E°56, = rie —Vvo, ) +O 5 = ety ag (3) 


where s*,* ,y;, — creep deformations; E — modulus of elasticity; v — Poisson's ratio; o;, oy, Try — values of stress 


J 


components in the corresponding directions. 
Having expressed the stress components in (3) through deformations, we write down the physical relations in the 


inverse form: 


o,= ; 
l-v 


E E 
(e, +Vvé, —(e% +ve',)); Oo, = al + VE, -(e% +ve));t,, = Tew -y%,). (4) 


The relationship of internal force factors and stresses is determined by integral relations: 


hl2 hl2 hl2 hi2 h/2 Al2 
N,= J 6,dz;N,= | o,dz;S= | t,dz;M,= | o,zdz;M,= | 0,2zdz;H = | t,,zdz, (5) 


—h/2 —h/2 —h/2 —h/2 —h/2 -h/2 


where N, and N, — linear longitudinal forces; S — linear shear forces, M, and M, — linear bending moments; 


H — linear torques; h — plate thickness. 
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Next, we substitute (1) in (4), as well as (4) in (5). As a result, we get: 


(6) 


where D = ea cylindrical rigidity of the plate, N* = z ‘a (e% + ve" )de 
12(1-v?) ? e Been l-v? -i2> * pes 
E ap E al2 
N= s\ + ve }dz, S* = ——~ * dz, 
*  J-y2 a . ) Mi+v) ft 
E Al2 E hi2 E h/2 
M*= e. + ve" )zdz, M* = e* + ve* )zdz, H* =—_——_ | _ y-,zdz 
*  J-y? Lf . :) *  J-v? LS :) ai+v) b 


Values N°, N*,S*,M*,M*, H* have the dimension of internal forces and determine the contribution of creep 


deformations to the redistribution of forces. 
Static equations of the flexible plate theory have the form [19]: 
ON, 
ON, ,_8 = 9-5 4 LG 
ox oy ox oy 
2 2 oM , Oo? 0? (wt 0? (wt 
pe ge yy ey CG) Mi) egg (wi) Ma 
Ox? Oxdy = Oy? : Ox? oy? Oxdy 
Here, gq — normal load on the surface of the plate, which is zero in this problem. 
It is possible to satisfy the first two static equations using the Airy stress functions, introduced by the formulas: 


2 2 2. 
ge oe Og r 
ae 


5) 


(7) 
q=0. 


Oxdy 


After substituting the last three equalities from (6) into the last static equation in (7) and taking into account (8), we 
obtain the first resolving equation: 
C20 02 (w+ Wy) C20 02 (w+ Wy) ) C20 0? (w+ Wy) 


DV‘w=q+q* + ————+ ———- 
one Ox? Oy? oy? Ox? Oxdy  OxOy ) 
2 * 2L7* 02M * 
where q* aM: oo 7 + |, 
Ox? Oxdy = Oy? 
To obtain the second resolving equation, it is required to express the deformation of the median surface from (6): 
N,-vN,+N*-vN- 
69 =— Sc Basia aint ee as jane vN* |; 
Eh Eh\ oy? Ox? ° 
N,-vN_+N*-—vN* 2 2 
g0 =—+__+*__*___* = OY a cage ane : (10) 
: Eh Eh\ ox? Oy? : 
2(1+ 2(1l+v 20 
go 2 ogg sO, | 
Eh Eh OxOy 
Substituting (10) into (2), we get: 
2 2 
| up = 6? (w+w) Ow, (wt my) (wt my) Omg OMG (tv) os + 
Eh Oxdy OxOy Ox? oy? ox? Oy? Eh Oxdy (11) 
O2N* O?N)\ @2Nt ea 
Vv x4 ; : 
Ox? oy? oy? Ox? 


Thus, for the problem under consideration, a system of resolving equations is obtained from two fourth-order 
differential equations (9) and (11). Equations (9) and (11) are nonlinear. In the resulting equations, values ® and w are 
functions of the coordinates x, y, and time ¢. Explicitly, there is no time in these equations, the time dependence is laid 
down in creep deformations ¢:, €;,, y%,, which are taken into account by the introduction of integral quantities 


Nt, Nt, S*, Mt, M*, He. 
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For the calculation scheme shown in Figure 1, the boundary conditions are written as: 


2 20 
atx=0,x =a: N, eo pS e 0;w=0;M, =0; 
“Oy? Oxdy ; 
oP Foti) (12) 
aty=0,y=b:N, = 0;S 0;w=0;M,, =0. 
Ox? Oxdy ; 


Equation (11) for small displacements, in the case of a plate made of elastic material, is a biharmonic equation that is 
used to solve the plane problem of the theory of elasticity in stresses. Frame analogy can serve as the boundary conditions 
for the stress function for the biharmonic equation. The plate contour is considered as a frame, and the stress function on 
the contour are equal to the bending moment / in it, and its derivative along the normal to the contour is the longitudinal 
force N. Plots M and N in the frame can be constructed in one of the basic systems of the force method (BSFM). The 
basic system, as well as the diagrams of the bending moment and longitudinal force in the frame are shown in Figure 2. 


Fig. 2. BSFM and diagrams of the bending moment and longitudinal forces 


If the vertical displacements do not exceed a quarter of the plate thickness, then it is possible to assume the forces in 
the middle surface independent of the coordinates x and y (N, =—p,N, =S =0) and use the linearized equation for 


calculations: 


(13) 


The analytical solution to the system of equations (9) and (11) is associated with great difficulties. The authors propose 
to solve this system numerically. The finite-difference method (FDM) was used in combination with the method of 
successive approximations. The Euler method was applied to determine creep deformation in the time domain. As the 
first stage, the solution for the elastic plate was performed. Load p was applied stepwise with small steps. At the initial 
load values, deflections w: were calculated by solving the simplified equation (13). Then, the calculated values w: were 
substituted into the differential equation (11). This provided determining the stress function. The next step was to solve 


the differential equation (9) using the known values of function ®, which made it possible to determine the nodal values 


of deflections w; . After that, the values (11 w, = (w, +wy ) /2 were substituted into equation (11). The iterative process 
was repeated at each step until the relative discrepancy with the norms of the vectors of the nodal values of deflections 
w, and w! was greater than the specified value (the authors assumed it to be equal to 0.1 %). For the second load step, 
the initial value w, in each node was the final result obtained in the first step. The calculation technique in the time 


domain, taking into account creep, was similar. The time interval, at which the process was investigated, was divided into 
steps At. In the case of setting the law of viscoelastic deformation in differential form, the values of creep deformations 


at step ¢ + At were calculated on the basis of the known rate of their growth at time ¢ using the Euler approximation: 


eo. =ert At. (14) 


The block diagram of the creep calculation algorithm is shown in Figure 3. 
Note that the system of equations (9) and (11) provides using schemes of a higher order of accuracy, e.g., the fourth- 


order Runge-Kutta method. At the same time, to achieve the same accuracy of the results, you can set noticeably large 
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time steps. However, with an increase in the step, there is a chance of not catching the effects of unsteady creep at the 
initial moments of time. And with the same time step, the Runge-Kutta method, in comparison to the Euler method, 


requires four times more operations to be performed. 


Solution to the elastic problem 
i=]t=0;e =e) = Vy =0;0=0,,, 


Output of results 


Fig. 3. Block diagram of the creep calculation algorithm 
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Research Results. A polyvinyl chloride polymer plate with dimensions a =2 m, b= 2 m, h=1 cm at E = 1480 MPa, 
v = 0.3 was considered. The nonlinear Maxwell-Gurevich equation was adopted as the law determining the rate of creep 


strain growth: 


0&;, _ iy a _—— 
2E=X,V, J=X); 
oto 
os 6, )- £€;5 
Si =5(5) -60 Wr 2b ip (15) 
* * is * 3 * 
aa] =No ol * ST iis = (o,, 0, ) E,£,, 2 
m 2 max 
where 6, — Kronecker symbol; 5,=J,/3, J,=0,+0, — first invariant of the stress tensor; nj, E. and 


m* — rheological parameters of the material, called the initial relaxation viscosity, the modulus of high elasticity, and the 


modulus of velocity. 
Indices rr in formula (15) indicate the direction of the primary stresses. 
For PVC, the authors took the values of rheological parameters from [20]: E.=5.99-10? MPa, m* = 12.6 MPa, 


ni, =5.44-10’ MPa:s. The shape of the initial deflection wo(x,v) was taken in accordance with the first form of the loss of 


stability of the plate made of elastic material: 
Wy (x9) =fysinsin =. (16) 
a 


For a plate made of elastic material without initial imperfections, the critical load, in the case of the whole ratio of 
sides a/b, was determined from the formula [19]: 
4n?D 
be 
To verify the developed calculation algorithm, the first step was to solve the elastic test problem and compare the 
results to the calculation in the finite element LIRA-CAD package (Fig. 4). The value of the arrow of the initial 
deflection fo was set to 0.15 mm. The grid size when using the FDM was 20x20, the number of load steps was 200. 
When calculating in the LIRA-CAD PC, the plate was divided by triangular finite elements with a triangulation step 
of 0.1 m. The load step size was assumed to be the same as when using the FDM. The value of the critical load for the 
elastic plate, calculated from formula (17), was 1340 N/m. Table 1 shows a comparison of vertical displacements in 
the center of the plate for different load values obtained by the author's method and using the finite element method 
(FEM). The deflections calculated using two alternative methods are quite close, except for the load of 1330 N/m. The 
deviation at this load value can be explained by the fact that when approaching the critical load, the movements rush to 


(17) 


Pr = 


infinity. 


—— oD —>*> EI 6@w6@6S Oa ee 
0 0.0365 0.0456 0.911 1.37 1.82 2.28 2.73 3.19 3.65 


Nonlinear loading 1 
Isofields of displacements by Z(G) 
Units — mm 


= 
YelsX 


Fig. 4. Isofield of vertical displacements in LIRA-CAD PC (p = 1330 N/m) 
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Table 1 
Comparison of calculation results using the author's method and with the help of FEM 
w:103, mm 
p, N/m 
LIRA-CAD the authors’ method 

133 16 16 

266 37 37 

399 63 63 

532 98 99 

665 146 148 

798 218 221 

931 336 342 
1,064 562 578 
1,197 1,163 1229 

1,330 3,646 4,229 


In [21], the possibility of transition from the solution of the elastic problem of calculating plates to the solution at the 
end of the creep process is shown. The value of the long-term critical load p.. can be obtained by replacing the cylindrical 
stiffness D of the elastic plate with the long-term cylindrical stiffness D., which is determined from the formula: 

ahs 


(18) 


1 1 vol 
where a =—+—;B=—+ : 
EE, E 2E, 

For viscoelastic rods and round plates, it has been previously established that in case p<p., the growth of 
displacements in time slows down, and the deflection arrow comes to a finite value. If p =p., deflections grow at a 
constant rate. At p > p., the rate of deflection growth increases. 

The authors also analyzed the nature of the growth of deflections over time for p < po, p = pw and p > p~ for different 
values of the maximum initial loss fo. The deflection curves over time in the center of the plate at p = 0.9 px, p = p» and 


p = 1.1 p»are shown respectively in Figures 5—7. 


w, mm 


1.4 ; T Tr T T T T 7 T 


— f=0.025 mm 


1.2 —— f=0.05 mm 


— f=0.1 mm 


1.0 j—— fo=0.15 mm 
0.8 
0.6 


0.4 


0.2 


Fig. 5. Deflection curves over time in the center of the plate at p = 0.9p« 
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—_ f=0.025 mm 
— f=0.05 mm 
—— fo=0.1 mm 


—— fo=0.15 mm 


0 2 4 6 8 10 12 14 16 18 20 ¢, hour 


— fo=0.15 mm 
—— fo=0.1 mm 
—— f=0.05 mm 
—— fo=0.025 mm 


0 1 2 3 4 5 6 a] 8 9 10 ¢, hour 


Fig. 7. Deflection curves over time in the center of the plate at p = 1.1 po 


Discussion and Conclusion. It can be seen from Fig. 5 that at p < px, the deflection arrow always comes to the final 
value, regardless of the values of the initial imperfections. At the same time, at p > p., the pattern of deflection growth 
obtained in [21] occurs only with small displacements. When deflections reach values exceeding about a quarter of the 
plate thickness, the rate of deformation growth begins to decrease even in the case of loads exceeding the long-term 
critical one. It should also be noted that there is a complete absence of a section with an increasing rate of displacement 
growth for plates with large initial curvatures. The identified effects can be explained by the redistribution of efforts 
N,, N,,S in the middle surface. 

Summarizing the above, we can conclude that the vertical movements of the plate pivotally supported along the 
contour, under the action of a compressive load on one axis, always come to a final value if the load does not exceed the 
instantaneous critical one. In other words, with the considered fastening and loading, the plate is in stable equilibrium 


under creep conditions. 


The obtained equations and the calculation algorithm make it possible to calculate plates made of arbitrary viscoelastic 


materials for any fastening options. The law of the relationship between stresses and creep deformations can also be set 


arbitrarily. 
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